
rm(list=ls())
setwd("/hpc/group/pbenfeylab/CheWei/")


suppressMessages(library(Seurat))
suppressMessages(library(tidyverse))
suppressMessages(library(dplyr))
suppressMessages(library(cowplot))
suppressMessages(library(data.table))
suppressMessages(library(RColorBrewer))
suppressMessages(library(OneR))

Warning message:
“package ‘ggplot2’ was built under R version 4.2.3”
Warning message:
“package ‘dplyr’ was built under R version 4.2.3”


sessionInfo() # New DCC env

R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Stream 8

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/ch416/miniconda3/envs/seu4/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=C.UTF-8           LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] OneR_2.2           RColorBrewer_1.1-3 data.table_1.14.8  cowplot_1.1.1     
 [5] forcats_0.5.2      stringr_1.5.0      dplyr_1.1.1        purrr_1.0.1       
 [9] readr_2.1.3        tidyr_1.3.0        tibble_3.2.1       ggplot2_3.4.2     
[13] tidyverse_1.3.2    SeuratObject_4.1.3 Seurat_4.1.1.9001 

loaded via a namespace (and not attached):
  [1] readxl_1.4.1           uuid_1.1-0             backports_1.4.1       
  [4] plyr_1.8.8             igraph_1.4.2           repr_1.1.4            
  [7] lazyeval_0.2.2         sp_1.6-0               splines_4.2.2         
 [10] RcppHNSW_0.4.1         listenv_0.9.0          scattermore_0.8       
 [13] digest_0.6.31          htmltools_0.5.5        fansi_1.0.4           
 [16] magrittr_2.0.3         tensor_1.5             googlesheets4_1.0.1   
 [19] cluster_2.1.4          ROCR_1.0-11            tzdb_0.3.0            
 [22] globals_0.16.2         modelr_0.1.10          matrixStats_0.63.0    
 [25] timechange_0.1.1       spatstat.sparse_3.0-1  colorspace_2.1-0      
 [28] rvest_1.0.3            ggrepel_0.9.3          haven_2.5.1           
 [31] crayon_1.5.2           jsonlite_1.8.4         progressr_0.13.0      
 [34] spatstat.data_3.0-1    survival_3.4-0         zoo_1.8-12            
 [37] glue_1.6.2             polyclip_1.10-4        gtable_0.3.3          
 [40] gargle_1.2.1           leiden_0.4.3           future.apply_1.10.0   
 [43] abind_1.4-5            scales_1.2.1           DBI_1.1.3             
 [46] spatstat.random_3.1-4  miniUI_0.1.1.1         Rcpp_1.0.10           
 [49] viridisLite_0.4.1      xtable_1.8-4           reticulate_1.28       
 [52] htmlwidgets_1.6.2      httr_1.4.5             ellipsis_0.3.2        
 [55] ica_1.0-3              pkgconfig_2.0.3        uwot_0.1.14           
 [58] dbplyr_2.2.1           deldir_1.0-6           utf8_1.2.3            
 [61] tidyselect_1.2.0       rlang_1.1.0            reshape2_1.4.4        
 [64] later_1.3.0            munsell_0.5.0          cellranger_1.1.0      
 [67] tools_4.2.2            cli_3.6.1              generics_0.1.3        
 [70] broom_1.0.2            ggridges_0.5.4         evaluate_0.20         
 [73] fastmap_1.1.1          goftest_1.2-3          fs_1.6.1              
 [76] fitdistrplus_1.1-8     RANN_2.6.1             pbapply_1.7-0         
 [79] future_1.32.0          nlme_3.1-162           mime_0.12             
 [82] xml2_1.3.3             compiler_4.2.2         plotly_4.10.1         
 [85] png_0.1-8              spatstat.utils_3.0-2   reprex_2.0.2          
 [88] stringi_1.7.12         RSpectra_0.16-1        lattice_0.21-8        
 [91] IRdisplay_1.1          Matrix_1.5-4           vctrs_0.6.2           
 [94] pillar_1.9.0           lifecycle_1.0.3        spatstat.geom_3.1-0   
 [97] lmtest_0.9-40          RcppAnnoy_0.0.20       irlba_2.3.5.1         
[100] httpuv_1.6.9           patchwork_1.1.2        R6_2.5.1              
[103] promises_1.2.0.1       KernSmooth_2.23-20     gridExtra_2.3         
[106] parallelly_1.35.0      codetools_0.2-19       fastDummies_1.6.3     
[109] MASS_7.3-58.3          assertthat_0.2.1       withr_2.5.0           
[112] sctransform_0.3.5      parallel_4.2.2         hms_1.1.2             
[115] grid_4.2.2             IRkernel_1.3.1.9000    googledrive_2.0.0     
[118] Rtsne_0.16             pbdZMQ_0.3-8           spatstat.explore_3.1-0
[121] shiny_1.7.4            lubridate_1.9.0        base64enc_0.1-3


pp.os.genes <- as.character(read.table("./CW_data/escoring/Root_sc/Rice_Protoplasting_DEgene_FC2_list.txt", header=F)$V1)
pp.os.genes <- gsub(pattern = "_", replacement = "-", x = pp.os.genes)
mt.os.genes <- as.character(read.table("./Meta_data/MSU7_mt_name.txt", header=F)$V1)
cp.os.genes <- as.character(read.table("./Meta_data/MSU7_cp_name.txt", header=F)$V1)
os2at <- read.csv("./CW_data/escoring/Root_sc/PLAZA_rice_arab_top_curated.csv", header=T)


use.sample <- c("tz2","tz1","sc_108","sc_109","sc_7","sc_115","sc_116","sc_192","sc_193","sc_194","sc_195","sc_196")
use.sample


read_seu <- function(dir,sample.name) { 
  #bscs <- read.csv("../../proj_sc/cbpsc/Benfey_single_cell-Samples.csv", na.strings=c("","NA"))
  #bscs <- bscs %>% select(c('sample','name','source','genotype','transgene','treatment','age','timepoint','rep','target_cells','date','seq_run')) %>% filter(sample==sample.name)
  #bscs$date <- gsub('^([0-9]{4})([0-9]{2})([0-9]+)$', '\\1-\\2-\\3', bscs$date)
  #bscs$target_cells <- prettyNum(bscs$target_cells, big.mark = ',')
  #bscs <- t(bscs)
  
  seu <- readRDS(dir) 
  #seu@misc$sample.meta.data <- list(bscs)
  seu@assays$spliced_RNA <- NULL
  seu@assays$spliced_SCT <- NULL
  seu@assays$unspliced_RNA <- NULL
  seu@assays$unspliced_SCT <- NULL
  return(seu)
}

rc.list <- list()

for (i in 1:length(use.sample))
{
  rc.list[[i]]<-read_seu(dir = paste0("./scRNA-seq/Seurat_Objects/",use.sample[i],"_COPILOT.rds"), sample.name = use.sample[i])
}

names(rc.list) <- use.sample

for (i in 1:length(use.sample))
{
  rc.list[[i]]<- suppressWarnings(suppressMessages(UpdateSeuratObject(rc.list[[i]])))
}


rc.features <- SelectIntegrationFeatures(object.list = rc.list, nfeatures = 40000)
length(rc.features)


rc.features <- rc.features[-match(pp.os.genes,rc.features)[!is.na(match(pp.os.genes,rc.features))]]
length(rc.features)


rc.features <- rc.features[-match(mt.os.genes,rc.features)[!is.na(match(mt.os.genes,rc.features))]]
length(rc.features)


rc.features <- rc.features[-match(cp.os.genes,rc.features)[!is.na(match(cp.os.genes,rc.features))]]
length(rc.features)


rc.list

$tz2
An object of class Seurat 
65686 features across 16754 samples within 2 assays 
Active assay: SCT (29172 features, 29172 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$tz1
An object of class Seurat 
56602 features across 15979 samples within 2 assays 
Active assay: SCT (24983 features, 24983 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_108
An object of class Seurat 
57342 features across 4995 samples within 2 assays 
Active assay: SCT (26163 features, 26163 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_109
An object of class Seurat 
58929 features across 4499 samples within 2 assays 
Active assay: SCT (26817 features, 26817 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_7
An object of class Seurat 
55557 features across 3042 samples within 2 assays 
Active assay: SCT (25311 features, 25311 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_115
An object of class Seurat 
59178 features across 3381 samples within 2 assays 
Active assay: SCT (26328 features, 26328 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_116
An object of class Seurat 
62181 features across 5188 samples within 2 assays 
Active assay: SCT (27719 features, 27719 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_192
An object of class Seurat 
60945 features across 6650 samples within 2 assays 
Active assay: SCT (27654 features, 27654 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_193
An object of class Seurat 
60565 features across 5296 samples within 2 assays 
Active assay: SCT (27308 features, 27308 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_194
An object of class Seurat 
60008 features across 4608 samples within 2 assays 
Active assay: SCT (27306 features, 27306 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_195
An object of class Seurat 
59194 features across 4802 samples within 2 assays 
Active assay: SCT (26979 features, 26979 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap

$sc_196
An object of class Seurat 
60276 features across 4549 samples within 2 assays 
Active assay: SCT (26304 features, 26304 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap


rc.list <- PrepSCTIntegration(object.list = rc.list, anchor.features = rc.features, verbose = TRUE)


rc.anchors <- suppressMessages(FindIntegrationAnchors(object.list = rc.list, normalization.method = "SCT", 
    anchor.features = rc.features, verbose = TRUE, reference=1))

Warning message in CheckDuplicateCellNames(object.list = object.list):
“Some cell names are duplicated across objects provided. Renaming to enforce unique cell names.”


# Running integration
rc.integrated <- IntegrateData(anchorset = rc.anchors, normalization.method = "SCT", verbose = TRUE)

Integrating dataset 2 with reference dataset

Finding integration vectors

Finding integration vector weights

Integrating data

Integrating dataset 3 with reference dataset

Finding integration vectors

Finding integration vector weights

Integrating data

Integrating dataset 4 with reference dataset

Finding integration vectors

Finding integration vector weights

Integrating data

Integrating dataset 5 with reference dataset

Finding integration vectors

Finding integration vector weights

Integrating data

Integrating dataset 6 with reference dataset

Finding integration vectors

Finding integration vector weights

Integrating data

Integrating dataset 7 with reference dataset

Finding integration vectors

Finding integration vector weights

Integrating data

Integrating dataset 8 with reference dataset

Finding integration vectors

Finding integration vector weights

Integrating data

Integrating dataset 9 with reference dataset

Finding integration vectors


# Run PCA
rc.integrated <- RunPCA(rc.integrated, npcs = 50, verbose = FALSE, approx = FALSE)

# Find nearest neighbors
rc.integrated <- FindNeighbors(rc.integrated, reduction = "pca",dims = 1:50)

# Clustering, notice that here we choose SLM algorithm with resoltion 5. Parameter "algorithm": 1 = original Louvain algorithm; 2 = Louvain algorithm with multilevel refinement; 3 = SLM algorithm; 4 = Leiden algorithm
rc.integrated <- FindClusters(rc.integrated, resolution = 0.5, algorithm = 3)


# Run 50D UMAP
rc.integrated <- RunUMAP(rc.integrated, reduction = "pca", dims = 1:50, n.components = 50)
rc.integrated@reductions$umap_50 <- rc.integrated@reductions$umap

# Run 3D UMAP
rc.integrated <- RunUMAP(rc.integrated, reduction = "pca", dims = 1:50, n.components = 3)
rc.integrated@reductions$umap_3D <- rc.integrated@reductions$umap

# Run 2D UMAP
rc.integrated <- RunUMAP(rc.integrated, reduction = "pca", dims = 1:50)
rc.integrated@reductions$umap_2D <- rc.integrated@reductions$umap


# Merge DoubletFinder results (optional) 
DF <- c()
for (i in grep("DF.classifications",colnames(rc.integrated@meta.data))){
    DF <- c(DF, rc.integrated@meta.data[,i][!is.na(rc.integrated@meta.data[,i])])
}
rc.integrated@meta.data$DF.classifications <- DF


rc.integrated[["umap"]]@cell.embeddings[,1] <- rc.integrated[["umap"]]@cell.embeddings[,1]*-1
#rc.integrated[["umap"]]@cell.embeddings[,2] <- rc.integrated[["umap"]]@cell.embeddings[,2]*-1
u2 <- rc.integrated@reductions$umap@cell.embeddings[,1]
u1 <- rc.integrated@reductions$umap@cell.embeddings[,2]
rc.integrated@reductions$umap@cell.embeddings[,1] <- u1
rc.integrated@reductions$umap@cell.embeddings[,2] <- u2


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, label=TRUE)+NoLegend()


options(repr.plot.width=10, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "orig.ident", cols= rainbow(12))


saveRDS(rc.integrated, file = "./scRNA-seq/Integrated_Objects/Rice_12S_WT_atlas_20221117_seu3.rds")


rc.integrated$log_nCount_RNA <- log(rc.integrated$nCount_RNA)


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "log_nCount_RNA")


load("./Meta_data/Root_bulk_rice_curated.RD")


nrow(time_MSU)


head(time_MSU)


# Extract matrix of integrated-batch-corrected expression value 
rc <- as.matrix(rc.integrated@assays$integrated@data)


# Merge the reference expression profile with the normalized expression matrix of our sample  

merge.rownames <- function (x,y){
  dat <- merge(x = x, y = y, by = "row.names")
  rownames(dat) <- dat$Row.names
  dat <- dat[,-1]
  return(dat)
}

time <- Reduce(merge.rownames, list(time_MSU,rc))


## Number of genes
nrow(time)


head(time)


# Prepare customized label name (optional)
time_label=c("Meristem","Elongation","Maturation1","Maturation2")


# Calculating the correlation coefficient of each cell to each reference expression profile and annotate the cell as the label that it has the highest correlation coefficient with.  
time_stat <- suppressWarnings(sapply(5:ncol(time), function(i) sapply(1:4, function(j) cor.test(time[,i],time[,j],method = "pearson")[c(3,4)])))
time_cor <- time_stat[seq(2,nrow(time_stat),2),]
time_pvalue <- time_stat[seq(1,nrow(time_stat)-1,2),]
time_max <- sapply(1:(ncol(time)-4), function(i) max(as.numeric(time_cor[,i])))
time_ident <- sapply(1:(ncol(time)-4), function(i) time_label[which(as.numeric(time_cor[,i])==max(as.numeric(time_cor[,i])))])
time_maxp <- sapply(1:(ncol(time)-4), function(i) as.numeric(time_pvalue[,i])[which(as.numeric(time_cor[,i])==max(as.numeric(time_cor[,i])))])
names(time_max) <- time_ident


# Store the annotation, correlation coefficient and the p-value in Seurat object
rc.integrated@meta.data$timezone.ID.P <- as.character(time_ident)
rc.integrated@meta.data$timezone.cor.P <- time_max
rc.integrated@meta.data$timezone.pvalue.P <- time_maxp

# In case there is cell with insufficient information for annotation, label them as "unknown"
rc.integrated@meta.data$timezone.ID.P[which(rc.integrated@meta.data$timezone.ID.P=='character(0)')]="unknown"


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "timezone.ID.P", order = c("Maturation2","Maturation1","Elongation","Meristem"),cols = c("#ebf8e3", "#51c8bd", "#009ac8", "#005fa8"))


table(rc.integrated$orig.ident)

sc_108 sc_109 sc_115 sc_116 sc_192 sc_193 sc_194 sc_195 sc_196   sc_7    tz1 
  4995   4499   3381   5188   6650   5296   4608   4802   4549   3042  15979 
   tz2 
 16754


rc.integrated$geno <- rc.integrated$orig.ident
rc.integrated$geno <- gsub("sc_7","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_108","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_109","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_115","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_116","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_192","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_193","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_194","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_195","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_196","X.Kitaake",rc.integrated$geno)
rc.integrated$geno <- gsub("tz1","ZH11",rc.integrated$geno)
rc.integrated$geno <- gsub("tz2","ZH11",rc.integrated$geno)


m1 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_7")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_7")])
m2 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_108")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_108")])
m3 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_109")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_109")])
m4 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_115")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_115")])
m5 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_116")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_116")])
m6 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_192")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_192")])
m7 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_193")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_193")])
m8 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_194")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_194")])
m9 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_195")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_195")])
m10 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_196")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="sc_196")])
m11 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="tz1")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="tz1")])
m12 <- rowSums(as.matrix(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="tz2")]))/ncol(rc.integrated@assays$SCT@data[,which(rc.integrated$orig.ident=="tz2")])

Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.2 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.1 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.3 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.6 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.3 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.1 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.2 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 1.1 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 3.9 GiB”
Warning message in asMethod(object):
“sparse->dense coercion: allocating vector of size 4.1 GiB”


dat <- data.frame(sc_7_Kitaake_WT_2days=m1, sc_108_Kitaake_WT_4days=m2, sc_109_Kitaake_WT_4days=m3, sc_115_Kitaake_WT_3days=m4, sc_116_Kitaake_WT_3days=m5,
                 sc_192_Kitaake_WT_4days=m6, sc_193_Kitaake_WT_4days=m7,sc_194_Kitaake_WT_4days=m8,sc_195_Kitaake_WT_4days=m9,sc_196_Kitaake_WT_4days=m10,
                 tz1_ZH11_WT_5days=m11,tz2_ZH11_WT_5days=m12)


options(repr.plot.width=12, repr.plot.height=12)
##The agglomeration method to be used. This should be (an unambiguous abbreviation of) one of "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).
##The distance measure to be used. This must be one of "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski". Any unambiguous substring can be given.
gplots::heatmap.2(cor(dat,method = "pearson"),col = brewer.pal(11,"YlOrRd"), 
          hclustfun = function(x) hclust(x, method = "median"), dendrogram="none",
          distfun = function(x) dist(x, method = "euclidean"), key.title = "Pearson correlation", trace = "none",cexCol = 1, cexRow = 1, margin=c(14,14),Colv=TRUE, Rowv = TRUE)

Warning message in brewer.pal(11, "YlOrRd"):
“n too large, allowed maximum for palette YlOrRd is 9
Returning the palette you asked for with that many colors
”


rc.integrated <- readRDS("./scRNA-seq/Integrated_Objects/Rice_12S_WT_atlas_20221117_seu3.rds")


head(rc.integrated@assays$integrated@data)


head(rc.integrated@assays$SCT@data)

  [[ suppressing 34 column names ‘AAACCCAAGAATACAC_1’, ‘AAACCCAAGAGTTCGG_1’, ‘AAACCCAAGGTAGCAC_1’ ... ]]

6 x 79743 sparse Matrix of class "dgCMatrix"
                                                                               
LOC-Os01g01010 .         . . . .         .         .        .         .        
LOC-Os01g01019 .         . . . .         .         .        .         .        
LOC-Os01g01030 .         . . . .         .         .        .         .        
LOC-Os01g01040 .         . . . 0.6931472 0.6931472 .        0.6931472 0.6931472
LOC-Os01g01050 .         . . . .         .         .        0.6931472 .        
LOC-Os01g01060 0.6931472 . . . .         .         1.791759 2.4849066 .        
                                                                      
LOC-Os01g01010 .         . .         .        . .        .         . .
LOC-Os01g01019 .         . .         .        . .        .         . .
LOC-Os01g01030 .         . .         .        . .        .         . .
LOC-Os01g01040 .         . 0.6931472 .        . 1.098612 0.6931472 . .
LOC-Os01g01050 .         . .         .        . .        .         . .
LOC-Os01g01060 0.6931472 . 0.6931472 2.079442 . .        .         . .
                                                                        
LOC-Os01g01010 .         . .        .        . . . .         . .        
LOC-Os01g01019 .         . .        .        . . . .         . .        
LOC-Os01g01030 .         . .        .        . . . .         . .        
LOC-Os01g01040 0.6931472 . .        .        . . . 0.6931472 . 0.6931472
LOC-Os01g01050 .         . .        .        . . . .         . .        
LOC-Os01g01060 .         . 1.386294 1.098612 . . . .         . .        
                                                                
LOC-Os01g01010 .         .         . . .        .         ......
LOC-Os01g01019 .         .         . . .        .         ......
LOC-Os01g01030 .         .         . . .        .         ......
LOC-Os01g01040 0.6931472 0.6931472 . . 1.098612 .         ......
LOC-Os01g01050 .         .         . . .        0.6931472 ......
LOC-Os01g01060 .         .         . . .        1.0986123 ......

 .....suppressing 79709 columns in show(); maybe adjust 'options(max.print= *, width = *)'
 ..............................


head(rc.integrated@assays$RNA@counts)

  [[ suppressing 34 column names ‘AAACCCAAGAATACAC_1’, ‘AAACCCAAGAGTTCGG_1’, ‘AAACCCAAGGTAGCAC_1’ ... ]]

6 x 79743 sparse Matrix of class "dgCMatrix"
                                                                               
LOC-Os01g01010 . . . . . .  .  . . . . . . . . . . . . .  . . . . . 1 . . . . .
LOC-Os01g01019 . . . . . .  .  . . . . . . . . . . . . .  . . . . . . . . . . .
LOC-Os01g01030 . . . . . .  .  . . . . . . . . . . . . .  . . . . . . . . . . .
LOC-Os01g01040 . . . 1 3 1  .  1 2 1 . 1 . . 1 2 . . 1 .  2 . . . . 4 1 2 1 1 .
LOC-Os01g01050 . . . . 1 .  .  1 . . . . . . . . . . . .  . . . . . . . . . . .
LOC-Os01g01060 1 . . . . . 13 15 . 2 . 1 3 . . . . . . . 18 1 . . . . . 1 . . .
                           
LOC-Os01g01010 . . . ......
LOC-Os01g01019 . . . ......
LOC-Os01g01030 . . . ......
LOC-Os01g01040 . 2 . ......
LOC-Os01g01050 . . 1 ......
LOC-Os01g01060 . . 2 ......

 .....suppressing 79709 columns in show(); maybe adjust 'options(max.print= *, width = *)'
 ..............................


## Normalized values (log normalization)
head(rc.integrated@assays$RNA@data)

  [[ suppressing 34 column names ‘AAACCCAAGAATACAC_1’, ‘AAACCCAAGAGTTCGG_1’, ‘AAACCCAAGGTAGCAC_1’ ... ]]

6 x 79743 sparse Matrix of class "dgCMatrix"
                                                                               
LOC-Os01g01010 . . . . . .  .  . . . . . . . . . . . . .  . . . . . 1 . . . . .
LOC-Os01g01019 . . . . . .  .  . . . . . . . . . . . . .  . . . . . . . . . . .
LOC-Os01g01030 . . . . . .  .  . . . . . . . . . . . . .  . . . . . . . . . . .
LOC-Os01g01040 . . . 1 3 1  .  1 2 1 . 1 . . 1 2 . . 1 .  2 . . . . 4 1 2 1 1 .
LOC-Os01g01050 . . . . 1 .  .  1 . . . . . . . . . . . .  . . . . . . . . . . .
LOC-Os01g01060 1 . . . . . 13 15 . 2 . 1 3 . . . . . . . 18 1 . . . . . 1 . . .
                           
LOC-Os01g01010 . . . ......
LOC-Os01g01019 . . . ......
LOC-Os01g01030 . . . ......
LOC-Os01g01040 . 2 . ......
LOC-Os01g01050 . . 1 ......
LOC-Os01g01060 . . 2 ......

 .....suppressing 79709 columns in show(); maybe adjust 'options(max.print= *, width = *)'
 ..............................


## How does a Seurat object look like?
rc.integrated

An object of class Seurat 
96918 features across 79743 samples within 3 assays 
Active assay: integrated (19395 features, 19395 variable features)
 2 other assays present: RNA, SCT
 5 dimensional reductions calculated: pca, umap, umap_50, umap_3D, umap_2D


## Annotations
head(rc.integrated@meta.data)


colnames(rc.integrated@meta.data)


table(rc.integrated$orig.ident)

sc_108 sc_109 sc_115 sc_116 sc_192 sc_193 sc_194 sc_195 sc_196   sc_7    tz1 
  4995   4499   3381   5188   6650   5296   4608   4802   4549   3042  15979 
   tz2 
 16754


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "orig.ident", cols= rainbow(12))


## Plot cell type
order <- c("Stem Cell Niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno , levels = order[sort(match(unique(rc.integrated$celltype.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols = color)+ggtitle("celltype.anno")


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno",order = c("Root Cap","Maturation","Elongation","Meristem"),cols = c("#e5f5e0", "#a1d99b", "#31a354", "#3182bd"))


zscore <- function(x){(x-mean(x))/sd(x)}
known.good.markers <- read.csv("./Meta_data/Potential_good_markers_for_annotation_20221129_Final.csv")
known.good.markers$Gene.ID <- gsub("^LOC_","LOC-",known.good.markers$Gene.ID)


known.good.markers


#known.good.markers$Cell.type <- gsub("Epidermis, Root cap", "Root Cap", known.good.markers$Cell.type, ignore.case = FALSE, perl = FALSE,fixed = T, useBytes = FALSE);
#known.good.markers$Celltype <- gsub("Epidermis", "Atrichoblast", known.good.markers$Celltype, ignore.case = FALSE, perl = FALSE,fixed = T, useBytes = FALSE);
#known.good.markers$Cell.type <- gsub("Early Metaxylem", "Xylem", known.good.markers$Cell.type, ignore.case = FALSE, perl = FALSE,fixed = T, useBytes = FALSE);
#known.good.markers$Celltype <- gsub("Late Metaxylem", "Xylem", known.good.markers$Celltype, ignore.case = FALSE, perl = FALSE,fixed = T, useBytes = FALSE);
known.good.markers$Cell.type <- gsub("Meristem", "Stem cell niche", known.good.markers$Cell.type, ignore.case = FALSE, perl = FALSE,fixed = T, useBytes = FALSE);
#known.good.markers$Cell.type <- gsub("QC", "Stem cell niche", known.good.markers$Cell.type, ignore.case = FALSE, perl = FALSE,fixed = T, useBytes = FALSE);
known.good.markers$Cell.type <- gsub("Root Cap junction", "Root Cap", known.good.markers$Cell.type, ignore.case = FALSE, perl = FALSE,fixed = T, useBytes = FALSE);


known.good.markers


## No markers not shared
which(is.na(match(known.good.markers$Gene.ID,rownames(rc.integrated@assays$integrated@data))))


#known.good.markers <- known.good.markers[-which(is.na(match(known.good.markers$Gene.ID,rownames(rc.integrated@assays$integrated@data)))),]


known.good.markers <- known.good.markers[-which(known.good.markers$Cell.type == "Xylem"),]


known.good.markers


as.character(unique(known.good.markers$Cell.type))


msc <- c()
for (i in as.character(unique(known.good.markers$Cell.type))){
    if (length(known.good.markers[which(known.good.markers$Cell.type== i),]$Gene.ID)>1){
    msc <- cbind(msc, as.numeric(apply(apply(rc.integrated@assays$integrated@data[known.good.markers[which(known.good.markers$Cell.type== i),]$Gene.ID,], 1, zscore), 1, mean)))       
    } else {
    msc <- cbind(msc, as.numeric(zscore(rc.integrated@assays$integrated@data[known.good.markers[which(known.good.markers$Cell.type== i),]$Gene.ID,])))      
    }

}
colnames(msc) <- as.character(unique(known.good.markers$Cell.type))
rownames(msc) <- colnames(rc.integrated)


# Crude
# Find clusters, here we choose Leiden clustering algorithm with resolution 2. Parameter "algorithm": 1 = original Louvain algorithm; 2 = Louvain algorithm with multilevel refinement; 3 = SLM algorithm; 4 = Leiden algorithm
DefaultAssay(rc.integrated) <- "integrated"
suppressMessages(suppressWarnings(
  rc.integrated <- FindClusters(rc.integrated, resolution = 2, algorithm = 3)
))

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 79743
Number of edges: 2436873

Running smart local moving algorithm...
Maximum modularity in 10 random starts: 0.8963
Number of communities: 62
Elapsed time: 72 seconds


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", label=TRUE)+NoLegend()


anno <- rc.integrated$seurat_clusters
for (i in unique(rc.integrated$seurat_clusters)){
    if (max(apply(msc[which(rc.integrated$seurat_clusters==i),],2,mean))>0){
        ct <- names(which.max(apply(msc[which(rc.integrated$seurat_clusters==i),],2,mean)))
    } else {
        ct <- "NA"
    }
        anno <- gsub(paste0("^",i,"$"), ct, anno, ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
}


rc.integrated$score.crude.anno <- anno


# Plot marker annotation
order <- c("Stem cell niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "NA")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$score.crude.anno <- factor(rc.integrated$score.crude.anno , levels = order[sort(match(unique(rc.integrated$score.crude.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$score.crude.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "score.crude.anno", cols = color)+ggtitle("Z-Score Annotation Crude")


# Find clusters, here we choose Leiden clustering algorithm with resolution 0.5. Parameter "algorithm": 1 = original Louvain algorithm; 2 = Louvain algorithm with multilevel refinement; 3 = SLM algorithm; 4 = Leiden algorithm
DefaultAssay(rc.integrated) <- "integrated"
suppressMessages(suppressWarnings(
  rc.integrated <- FindClusters(rc.integrated, resolution = 300, algorithm = 3)
))

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 79743
Number of edges: 2436873

Running smart local moving algorithm...
Maximum modularity in 10 random starts: 0.4437
Number of communities: 1929
Elapsed time: 49 seconds


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap")+NoLegend()


anno <- rc.integrated$seurat_clusters
for (i in unique(rc.integrated$seurat_clusters)){
    if (max(apply(msc[which(rc.integrated$seurat_clusters==i),],2,mean))>0){
        ct <- names(which.max(apply(msc[which(rc.integrated$seurat_clusters==i),],2,mean)))
    } else {
        ct <- "NA"
    }
        anno <- gsub(paste0("^",i,"$"), ct, anno, ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
}


rc.integrated$score.anno <- anno


# Plot marker annotation
order <- c("Stem cell niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "NA")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$score.anno <- factor(rc.integrated$score.anno , levels = order[sort(match(unique(rc.integrated$score.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$score.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "score.anno", cols = color)+ggtitle("Z-Score Annotation")


## QC marker
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os01g63510")


## Pericycle marker
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os10g39980")


## New Root Cap Marker
DefaultAssay(rc.integrated) <- "SCT"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os03g18130")


## New Root Cap Marker 
DefaultAssay(rc.integrated) <- "integrated"
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os03g18130")


## Xylem markers (Marker from Zhang et al.)
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os08g38170")


## Xylem markers (Marker from Benfey lab)
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os01g48130")


## Xylem markers (Marker from Benfey lab)
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os09g25490")


## Xylem markers (Marker from Benfey lab)
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os01g54620")


## Xylem markers (Marker from Benfey lab)
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os10g32980")


## New trichoblast markers (Marker from Zhang)
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os10g42750")


## New trichoblast markers (Marker from Ding)
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os06g08500")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os01g50820")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os01g64840")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os03g19990")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os01g03530")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os01g11750")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os03g04210")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os04g45290")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os05g38770")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os06g48050")


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os10g03400")


## Root Cap juctions
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, feature= "LOC-Os10g31660")

Warning message:
“Found the following features in more than one assay, excluding the default. We will not include these in the final dataframe: LOC-Os10g31660”
Warning message in FetchData(object = object, vars = c(dims, "ident", features), :
“The following requested variables were not found: LOC-Os10g31660”

Error: None of the requested features were found: LOC-Os10g31660 in slot data
Traceback:

1. FeaturePlot(rc.integrated, feature = "LOC-Os10g31660")
2. stop("None of the requested features were found: ", paste(features, 
 .     collapse = ", "), " in slot ", slot, call. = FALSE)


dat <- data.frame(rc.integrated$score.anno, rc.integrated$score.crude.anno, rc.integrated$timezone.ID.P)


rc.integrated$consensus.anno <- apply(dat,1,function(x){if(x[1]=="NA" & x[2]=="NA"){"Unknown"}else if(x[1]==x[2]){x[1]}else if(x[1]=="Trichoblast" & x[2]=="Atrichoblast"){"Trichoblast"}else if(x[1]=="Vascular Tissue" & x[2]=="NA"){"Vascular Tissue"}
    else if(x[1]=="Metaxylem" & x[2]=="Phloem"){"Metaxylem"}else if(x[1]=="Cortex" & x[2]=="Sclerenchyma"){"Cortex"}
        else if(x[1]=="Exodermis" & x[2]=="Sclerenchyma"){"Exodermis"}else if(x[1]=="Exodermis" & x[2]=="Endodermis"){"Exodermis"}
            else if(x[1]=="Cortex" & x[2]=="Endodermis"){"Cortex"}else if(x[1]=="Pericycle" & x[2]=="Endodermis"){"Pericycle"}else if(x[1]=="Phloem" & x[2]=="Endodermis"){"Phloem"}
                    else if(x[1]=="Metaxylem" & x[2]=="Endodermis"){"Metaxylem"}else if (x[1]=="NA"|x[2]=="NA"){"Unknown"}else if(x[3]=="Maturation1"|x[3]=="Maturation2"){x[2]}else {"Unknown"}})


order <- c("Stem cell niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$consensus.anno <- factor(rc.integrated$consensus.anno , levels = order[sort(match(unique(rc.integrated$consensus.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$consensus.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "consensus.anno", cols = color)+ggtitle("consensus.anno")


rc.integrated$new_ref <- paste0(gsub("Maturation.$","Maturation",rc.integrated$timezone.ID.P), "_",rc.integrated$consensus.anno)


table(rc.integrated$new_ref)

   Elongation_Atrichoblast          Elongation_Cortex 
                      1091                       8944 
     Elongation_Endodermis       Elongation_Exodermis 
                        43                       4218 
      Elongation_Metaxylem       Elongation_Pericycle 
                       104                          9 
         Elongation_Phloem        Elongation_Root Cap 
                        44                         20 
   Elongation_Sclerenchyma     Elongation_Trichoblast 
                      4716                         80 
        Elongation_Unknown    Maturation_Atrichoblast 
                      1275                       1467 
         Maturation_Cortex      Maturation_Endodermis 
                      1183                        302 
      Maturation_Exodermis       Maturation_Metaxylem 
                      4792                       1545 
      Maturation_Pericycle          Maturation_Phloem 
                      1186                        338 
     Maturation_Protoxylem        Maturation_Root Cap 
                      3669                       4866 
   Maturation_Sclerenchyma Maturation_Stem cell niche 
                      1882                          1 
    Maturation_Trichoblast         Maturation_Unknown 
                      5007                       4586 
     Meristem_Atrichoblast            Meristem_Cortex 
                      1273                       4480 
       Meristem_Endodermis         Meristem_Exodermis 
                      1974                       1510 
        Meristem_Metaxylem         Meristem_Pericycle 
                       474                       4550 
           Meristem_Phloem          Meristem_Root Cap 
                      1078                        312 
     Meristem_Sclerenchyma   Meristem_Stem cell niche 
                      1299                       3021 
      Meristem_Trichoblast           Meristem_Unknown 
                         8                       8396


rc.integrated$new_ref <- gsub(".*Unknown$","Unknown",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub(".*Root Cap$","Root Cap",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub("Elongation_Stem cell niche","Unknown",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub("Maturation_Stem cell niche","Unknown",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub("^Meristem.*","Meristem",rc.integrated$new_ref)


table(rc.integrated$new_ref)

Elongation_Atrichoblast       Elongation_Cortex   Elongation_Endodermis 
                   1091                    8944                      43 
   Elongation_Exodermis    Elongation_Metaxylem    Elongation_Pericycle 
                   4218                     104                       9 
      Elongation_Phloem Elongation_Sclerenchyma  Elongation_Trichoblast 
                     44                    4716                      80 
Maturation_Atrichoblast       Maturation_Cortex   Maturation_Endodermis 
                   1467                    1183                     302 
   Maturation_Exodermis    Maturation_Metaxylem    Maturation_Pericycle 
                   4792                    1545                    1186 
      Maturation_Phloem   Maturation_Protoxylem Maturation_Sclerenchyma 
                    338                    3669                    1882 
 Maturation_Trichoblast                Meristem                Root Cap 
                   5007                   19667                    5198 
                Unknown 
                  14258


length(sort(table(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")])))


nnc <- length(sort(table(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")]))) 
# Extract integrated (batch-corrected) expression matrix
afm <- rc.integrated@assays$integrated@data

# Pool (average) expression values of each grouping
new_ref <- matrix(nrow=nrow(afm), ncol = 0)
for (i in 1:nnc) {
m <- afm[,which(rc.integrated$new_ref==as.character(unique(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")]))[i])]
new_ref <- cbind(new_ref, rowSums(m)/ncol(m))
}
colnames(new_ref) <- as.character(unique(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")]))
gene.var <- apply(new_ref,1,var)

# Select top 200 highly variable genes 
new_ref_sub <- new_ref[names(sort(gene.var,decreasing = TRUE)[1:200]),]

# Merge newly-built reference with atlas

merge.rownames <- function (x,y){
  dat <- merge(x = x, y = y, by = "row.names")
  rownames(dat) <- dat$Row.names
  dat <- dat[,-1]
  return(dat)
}

nr <- Reduce(merge.rownames, list(new_ref_sub,afm)) 
nr <- as.matrix(nr)


colnames(nr)[1:nnc+1]


# Correlation-based annotation using newly-built references
nr_label=colnames(new_ref)
nr_stat <- suppressWarnings(sapply((nnc+1):ncol(nr), function(i) sapply(1:nnc, function(j) cor.test(nr[,i],nr[,j],method = "pearson")[c(3,4)])))
nr_cor <- nr_stat[seq(2,nrow(nr_stat),2),]
nr_pvalue <- nr_stat[seq(1,nrow(nr_stat)-1,2),]                                                                                                        
nr_max <- sapply(1:(ncol(nr)-nnc), function(i) max(as.numeric(nr_cor[,i])))
nr_ident <- sapply(1:(ncol(nr)-nnc), function(i) nr_label[which(as.numeric(nr_cor[,i])==max(as.numeric(nr_cor[,i])))])
nr_maxp <- sapply(1:(ncol(nr)-nnc), function(i) as.numeric(nr_pvalue[,i])[which(as.numeric(nr_cor[,i])==max(as.numeric(nr_cor[,i])))])                     
names(nr_max) <- nr_ident
rc.integrated$final.anno.P <- as.character(nr_ident)
rc.integrated$final.anno.cor.P <- as.numeric(nr_max)
rc.integrated$final.anno.pvalue.P <- as.numeric(nr_maxp)


rc.integrated$celltype.anno <-  gsub("^Meristem_","",gsub("^Elongation_","",gsub("^Maturation_","",as.character(rc.integrated$final.anno.P))))


order <- c("Meristem", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno , levels = order[sort(match(unique(rc.integrated$celltype.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols = color)+ggtitle("celltype.anno")


rc.integrated$time.anno <-  gsub("_.*$","",as.character(rc.integrated$final.anno.P))


table(rc.integrated$time.anno)

Elongation Maturation   Meristem   Root Cap 
     36589      34421       5009       3724


options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno",order = c("Root Cap","Maturation","Elongation","Meristem"),cols = c("#e5f5e0", "#a1d99b", "#31a354", "#3182bd"))


options(repr.plot.width=10, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "geno")


rc.integrated$new_ref <- paste0(rc.integrated$time.anno, "_",rc.integrated$consensus.anno)
rc.integrated$new_ref <- gsub(".*Unknown$","Unknown",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub(".*Root Cap$","Root Cap",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub("^Root Cap_.*","Root Cap",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub("Elongation_Stem cell niche","Unknown",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub("Maturation_Stem cell niche","Unknown",rc.integrated$new_ref)
rc.integrated$new_ref <- gsub("^Meristem.*","Meristem",rc.integrated$new_ref)


table(rc.integrated$new_ref)

Elongation_Atrichoblast       Elongation_Cortex   Elongation_Endodermis 
                   2550                   12273                    1832 
   Elongation_Exodermis    Elongation_Metaxylem    Elongation_Pericycle 
                   5790                     500                    1271 
      Elongation_Phloem Elongation_Sclerenchyma  Elongation_Trichoblast 
                    675                    5404                     229 
Maturation_Atrichoblast       Maturation_Cortex   Maturation_Endodermis 
                   1277                    2233                     406 
   Maturation_Exodermis    Maturation_Metaxylem    Maturation_Pericycle 
                   4324                    1360                    4240 
      Maturation_Phloem   Maturation_Protoxylem Maturation_Sclerenchyma 
                    783                    3598                    2178 
 Maturation_Trichoblast                Meristem                Root Cap 
                   4865                    3122                    5737 
                Unknown 
                  15096


#rc.integrated$new_ref <- gsub("Meristem_Vascular Tissue","Unknown",rc.integrated$new_ref)


table(rc.integrated$new_ref)

Elongation_Atrichoblast       Elongation_Cortex   Elongation_Endodermis 
                   2550                   12273                    1832 
   Elongation_Exodermis    Elongation_Metaxylem    Elongation_Pericycle 
                   5790                     500                    1271 
      Elongation_Phloem Elongation_Sclerenchyma  Elongation_Trichoblast 
                    675                    5404                     229 
Maturation_Atrichoblast       Maturation_Cortex   Maturation_Endodermis 
                   1277                    2233                     406 
   Maturation_Exodermis    Maturation_Metaxylem    Maturation_Pericycle 
                   4324                    1360                    4240 
      Maturation_Phloem   Maturation_Protoxylem Maturation_Sclerenchyma 
                    783                    3598                    2178 
 Maturation_Trichoblast                Meristem                Root Cap 
                   4865                    3122                    5737 
                Unknown 
                  15096


length(sort(table(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")])))


nnc <- length(sort(table(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")])))


# Extract integrated (batch-corrected) expression matrix
afm <- rc.integrated@assays$integrated@data

# Pool (average) expression values of each grouping
new_ref <- matrix(nrow=nrow(afm), ncol = 0)
for (i in 1:nnc) {
m <- afm[,which(rc.integrated$new_ref==as.character(unique(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")]))[i])]
new_ref <- cbind(new_ref, rowSums(m)/ncol(m))
}
colnames(new_ref) <- as.character(unique(rc.integrated$new_ref[which(rc.integrated$new_ref!="Unknown")]))
gene.var <- apply(new_ref,1,var)

# Select top 200 highly variable genes 
new_ref_sub <- new_ref[names(sort(gene.var,decreasing = TRUE)[1:200]),]

# Merge newly-built reference with atlas

merge.rownames <- function (x,y){
  dat <- merge(x = x, y = y, by = "row.names")
  rownames(dat) <- dat$Row.names
  dat <- dat[,-1]
  return(dat)
}

nr <- Reduce(merge.rownames, list(new_ref_sub,afm)) 
nr <- as.matrix(nr)


colnames(nr)[1:nnc+1]


# Correlation-based annotation using newly-built references
nr_label=colnames(new_ref)
nr_stat <- suppressWarnings(sapply((nnc+1):ncol(nr), function(i) sapply(1:nnc, function(j) cor.test(nr[,i],nr[,j],method = "pearson")[c(3,4)])))
nr_cor <- nr_stat[seq(2,nrow(nr_stat),2),]
nr_pvalue <- nr_stat[seq(1,nrow(nr_stat)-1,2),]                                                                                                        
nr_max <- sapply(1:(ncol(nr)-nnc), function(i) max(as.numeric(nr_cor[,i])))
nr_ident <- sapply(1:(ncol(nr)-nnc), function(i) nr_label[which(as.numeric(nr_cor[,i])==max(as.numeric(nr_cor[,i])))])
nr_maxp <- sapply(1:(ncol(nr)-nnc), function(i) as.numeric(nr_pvalue[,i])[which(as.numeric(nr_cor[,i])==max(as.numeric(nr_cor[,i])))])                     
names(nr_max) <- nr_ident
rc.integrated$final.anno.P <- as.character(nr_ident)
rc.integrated$final.anno.cor.P <- as.numeric(nr_max)
rc.integrated$final.anno.pvalue.P <- as.numeric(nr_maxp)


rc.integrated$celltype.fine.anno <-  gsub("^Meristem_","",gsub("^Elongation_","",gsub("^Maturation_","",as.character(rc.integrated$final.anno.P))))
rc.integrated$celltype.fine.anno <- gsub("Meristem","Stem Cell Niche",rc.integrated$celltype.fine.anno)


order <- c("Stem Cell Niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.fine.anno <- factor(rc.integrated$celltype.fine.anno , levels = order[sort(match(unique(rc.integrated$celltype.fine.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.fine.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.fine.anno", cols = color)+ggtitle("celltype.fine.anno")


#rc.integrated$celltype.anno <-  gsub("Late Metaxylem","Xylem",gsub("Early Metaxylem","Xylem",gsub("^Meristem_","",gsub("^Elongation_","",gsub("^Maturation_","",as.character(rc.integrated$final.anno.P))))))
rc.integrated$celltype.anno <-  gsub("^Meristem_","",gsub("^Elongation_","",gsub("^Maturation_","",as.character(rc.integrated$final.anno.P))))
rc.integrated$celltype.anno <- gsub("Meristem","Stem Cell Niche",rc.integrated$celltype.anno)


order <- c("Stem Cell Niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno , levels = order[sort(match(unique(rc.integrated$celltype.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols = color)+ggtitle("celltype.anno")


rc.integrated$time.anno <-  gsub("_.*$","",as.character(rc.integrated$final.anno.P))


table(rc.integrated$time.anno)

Elongation Maturation   Meristem   Root Cap 
     34164      33967       7032       4580


options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno",order = c("Root Cap","Maturation","Elongation","Meristem"),cols = c("#e5f5e0", "#a1d99b", "#31a354", "#3182bd"))


saveRDS(rc.integrated, file = "./scRNA-seq/Integrated_Objects/Rice_12S_WT_atlas_20221117_seu3.rds")


## Celltype
cat("X.Kitaake")
table(rc.integrated$celltype.anno[which(rc.integrated$geno=="X.Kitaake")])/length(which(rc.integrated$geno=="X.Kitaake"))*100
cat("ZH11")
table(rc.integrated$celltype.anno[which(rc.integrated$geno=="ZH11")])/length(which(rc.integrated$geno=="ZH11"))*100
cat("sc_7")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_7")])/length(which(rc.integrated$orig.ident=="sc_7"))*100
cat("sc_108")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_108")])/length(which(rc.integrated$orig.ident=="sc_108"))*100
cat("sc_109")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_109")])/length(which(rc.integrated$orig.ident=="sc_109"))*100
cat("sc_115")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_115")])/length(which(rc.integrated$orig.ident=="sc_115"))*100
cat("sc_116")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_116")])/length(which(rc.integrated$orig.ident=="sc_116"))*100
cat("sc_192")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_192")])/length(which(rc.integrated$orig.ident=="sc_192"))*100
cat("sc_193")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_193")])/length(which(rc.integrated$orig.ident=="sc_193"))*100
cat("sc_194")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_194")])/length(which(rc.integrated$orig.ident=="sc_194"))*100
cat("sc_195")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_195")])/length(which(rc.integrated$orig.ident=="sc_195"))*100
cat("sc_196")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="sc_196")])/length(which(rc.integrated$orig.ident=="sc_196"))*100
cat("tz1")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="tz1")])/length(which(rc.integrated$orig.ident=="tz1"))*100
cat("tz2")
table(rc.integrated$celltype.anno[which(rc.integrated$orig.ident=="tz2")])/length(which(rc.integrated$orig.ident=="tz2"))*100

X.Kitaake

Stem Cell Niche    Sclerenchyma        Root Cap    Atrichoblast     Trichoblast 
       7.945118        8.951287        7.945118        3.839609        9.193789 
         Cortex      Endodermis          Phloem       Pericycle       Exodermis 
      21.214635       11.523080        2.646245        9.089555       10.321208 
     Protoxylem       Metaxylem 
       5.479685        1.850670

ZH11

Stem Cell Niche    Sclerenchyma        Root Cap    Atrichoblast     Trichoblast 
      10.072404       14.337213        2.581493        9.608041        3.528549 
         Cortex      Endodermis          Phloem       Pericycle       Exodermis 
      19.759875        8.309657        1.906333        5.899245       17.312804 
     Protoxylem       Metaxylem 
       4.332020        2.352366

sc_7

Stem Cell Niche    Sclerenchyma        Root Cap    Atrichoblast     Trichoblast 
       7.199211        6.344510        4.174885        2.103879        2.432610 
         Cortex      Endodermis          Phloem       Pericycle       Exodermis 
      28.172255       15.811966        2.235371       16.699540        6.278764 
     Protoxylem       Metaxylem 
       6.147272        2.399737

sc_108

Stem Cell Niche    Sclerenchyma        Root Cap    Atrichoblast     Trichoblast 
      10.090090        4.224224        7.027027        6.206206        5.325325 
         Cortex      Endodermis          Phloem       Pericycle       Exodermis 
      18.358358       13.093093        3.443443       17.337337        7.887888 
     Protoxylem       Metaxylem 
       2.562563        4.444444

sc_109

Stem Cell Niche    Sclerenchyma        Root Cap    Atrichoblast     Trichoblast 
       7.068237        4.378751        4.645477        3.734163       17.626139 
         Cortex      Endodermis          Phloem       Pericycle       Exodermis 
      16.203601       17.826184        2.089353       10.846855        9.402089 
     Protoxylem       Metaxylem 
       3.689709        2.489442


## Timezone
cat("X.Kitaake")
table(rc.integrated$time.anno[which(rc.integrated$geno=="X.Kitaake")])/length(which(rc.integrated$geno=="X.Kitaake"))*100
cat("ZH11")
table(rc.integrated$time.anno[which(rc.integrated$geno=="ZH11")])/length(which(rc.integrated$geno=="ZH11"))*100
cat("sc_7")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_7")])/length(which(rc.integrated$orig.ident=="sc_7"))*100
cat("sc_108")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_108")])/length(which(rc.integrated$orig.ident=="sc_108"))*100
cat("sc_109")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_109")])/length(which(rc.integrated$orig.ident=="sc_109"))*100
cat("sc_115")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_115")])/length(which(rc.integrated$orig.ident=="sc_115"))*100
cat("sc_116")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_116")])/length(which(rc.integrated$orig.ident=="sc_116"))*100
cat("sc_192")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_192")])/length(which(rc.integrated$orig.ident=="sc_192"))*100
cat("sc_193")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_193")])/length(which(rc.integrated$orig.ident=="sc_193"))*100
cat("sc_194")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_194")])/length(which(rc.integrated$orig.ident=="sc_194"))*100
cat("sc_195")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_195")])/length(which(rc.integrated$orig.ident=="sc_195"))*100
cat("sc_196")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="sc_196")])/length(which(rc.integrated$orig.ident=="sc_196"))*100
cat("tz1")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="tz1")])/length(which(rc.integrated$orig.ident=="tz1"))*100
cat("tz2")
table(rc.integrated$time.anno[which(rc.integrated$orig.ident=="tz2")])/length(which(rc.integrated$orig.ident=="tz2"))*100

X.Kitaake

Elongation Maturation   Meristem   Root Cap 
 38.162093  45.947671   7.945118   7.945118

ZH11

Elongation Maturation   Meristem   Root Cap 
 49.564660  37.781444  10.072404   2.581493

sc_7

Elongation Maturation   Meristem   Root Cap 
 36.653517  51.972387   7.199211   4.174885

sc_108

Elongation Maturation   Meristem   Root Cap 
 40.360360  42.522523  10.090090   7.027027

sc_109

Elongation Maturation   Meristem   Root Cap 
 39.430985  48.855301   7.068237   4.645477


pseudotime <- read.csv("./rice_ct_psudotime_SCT.csv")


rc.integrated$ct_pseudotime <- pseudotime$ct_pseudotime


options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(rc.integrated, features = "ct_pseudotime", pt.size=0.5)+ scale_colour_gradientn(colours = rev(brewer.pal(11,"Spectral")))

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.


df <- data.frame(celltype.anno = rc.integrated$celltype.anno, ct_pseudotime=rc.integrated$ct_pseudotime, pseudotime.group=rep("Unknown", ncol(rc.integrated)))


for (i in unique(df$celltype.anno)){
    if (i == "Stem Cell Niche") {
        df$pseudotime.group[which(df$celltype.anno==i)] <- "T0"
    } else {
        df$pseudotime.group[which(df$celltype.anno==i)] <- as.character(bin(df$ct_pseudotime[which(df$celltype.anno==i)], nbins=9, labels = paste0("T",seq(1,9,1)), method = "content")) 
    }
}


rc.integrated$pseudotime.group <- df$pseudotime.group


# Plot pseudotime group
nb.cols <- 10
mycolors <- colorRampPalette(brewer.pal(11, "Spectral"))(nb.cols)
order <- paste0("T",seq(0,9,1))
rc.integrated$pseudotime.group <- factor(rc.integrated$pseudotime.group, levels = order[sort(match(unique(rc.integrated$pseudotime.group),order))]) 
color <- mycolors[sort(match(unique(rc.integrated$pseudotime.group),order))]
options(repr.plot.width=10, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "pseudotime.group", pt.size=0.5, cols=color)


suppressMessages(library(plotly))


options(repr.plot.width=8, repr.plot.height=8)
rc.integrated@reductions$umap_2D <- rc.integrated@reductions$umap
rc.integrated@reductions$umap <- rc.integrated@reductions$umap_3D
order <- c("Stem Cell Niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno , levels = order[sort(match(unique(rc.integrated$celltype.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols = color)+ggtitle("celltype.anno")


u3 <- rc.integrated@reductions$umap@cell.embeddings[,1]
u1 <- rc.integrated@reductions$umap@cell.embeddings[,3]
rc.integrated@reductions$umap@cell.embeddings[,1] <- u1
rc.integrated@reductions$umap@cell.embeddings[,3] <- u3


int.name <- "Rice_Atlas"


rc.integrated@reductions$umap@cell.embeddings[,3] <- -rc.integrated@reductions$umap@cell.embeddings[,3]


# Make sure that the color match to the label it is supposed to color 
umap_1 <- rc.integrated[["umap"]]@cell.embeddings[,1]
umap_2 <- rc.integrated[["umap"]]@cell.embeddings[,2]
umap_3 <- rc.integrated[["umap"]]@cell.embeddings[,3]


# Prepare data for plotting
plot.data <- FetchData(object = rc.integrated, vars = c("UMAP_1", "UMAP_2", "UMAP_3", "celltype.anno", "time.anno", "orig.ident"))

# Make a column of row name identities (these will be your cell/barcode names)
plot.data$label <- paste(plot.data$orig.ident)

order <- c("Stem Cell Niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
p <- plot_ly(data = plot.data, 
        x = ~UMAP_1, y = ~UMAP_2, z = ~UMAP_3, 
        color = ~celltype.anno, 
        colors = color,
        type = "scatter3d", 
        mode = "markers", 
        marker = list(size = 3, width=2), # controls size of points
        text=~label, #This is that extra column we made earlier for which we will use for cell ID
        hoverinfo="text") #When you visualize your plotly object, hovering your mouse pointer over a point shows cell names
htmlwidgets::saveWidget(p, paste0(int.name,"_celltype.anno.html"))

order <- c("Maturation","Elongation","Meristem","Root Cap")
palette <- c("#31a354", "#a1d99b", "#e5f5e0", "#3182bd")
rc.integrated$time.anno <- factor(rc.integrated$time.anno, levels = order[sort(match(unique(rc.integrated$time.anno),order))])
color <- palette[sort(match(unique(rc.integrated$time.anno),order))]
p <- plot_ly(data = plot.data, 
        x = ~UMAP_1, y = ~UMAP_2, z = ~UMAP_3, 
        color = ~time.anno, 
        colors = color,
        type = "scatter3d", 
        mode = "markers", 
        marker = list(size = 3, width=2), # controls size of points
        text=~label, #This is that extra column we made earlier for which we will use for cell ID
        hoverinfo="text") #When you visualize your plotly object, hovering your mouse pointer over a point shows cell names
htmlwidgets::saveWidget(p, paste0(int.name,"_time.anno.html"))

Warning message in `[[.DimReduc`(object[[x]], cells, vars.use, drop = FALSE):
“The following embeddings are not present: UMAP_3”


options(repr.plot.width=8, repr.plot.height=8)
order <- c("Stem Cell Niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno , levels = order[sort(match(unique(rc.integrated$celltype.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols = color)+ggtitle("celltype.anno")


rc.integrated@reductions$umap <- rc.integrated@reductions$umap_2D


options(repr.plot.width=8, repr.plot.height=8)
order <- c("Stem Cell Niche", "Sclerenchyma", "Columella", "Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem", "Xylem", "Vascular Tissue","Pericycle","Exodermis","Epidermis", "Early Metaxylem", "Late Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#dd77ec", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno , levels = order[sort(match(unique(rc.integrated$celltype.anno),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols = color)+ggtitle("celltype.anno")


saveRDS(rc.integrated, file = "./Rice_Atlas_12WT_ref_tz2_seu3.rds")


rc.integrated <- readRDS("./Rice_Atlas_12WT_ref_tz2_seu3.rds")


rc.integrated <- UpdateSeuratObject(rc.integrated)

Validating object structure

Updating object slots

Ensuring keys are in the proper strucutre

Ensuring feature names don't have underscores or pipes

Object representation is consistent with the most current Seurat version


DefaultAssay(rc.integrated) <- "integrated"


rc.integrated

An object of class Seurat 
88280 features across 44144 samples within 3 assays 
Active assay: integrated (16838 features, 16838 variable features)
 2 other assays present: RNA, SCT
 5 dimensional reductions calculated: pca, umap, umap_50, umap_3D, umap_2D


# Need to change the directory accordingly 
label_transfer_init <- function(sample.name){
    seu <-read_seu(dir = paste0("../COPILOT_RDS_LENIENT_DR/",sample.name,"_COPILOT.rds"))
    lt.anchors <- FindTransferAnchors(reference = rc.integrated, query = seu, normalization.method = "SCT", npcs = 50 ,dims = 1:50, features = intersect(rownames(rc.integrated),rownames(seu)))
    predictions <- TransferData(anchorset = lt.anchors, refdata = rc.integrated$celltype.anno, dims = 1:50, weight.reduction = "pcaproject")
    seu <- AddMetaData(seu, metadata = predictions)
    seu@meta.data$celltype.anno <- seu@meta.data$predicted.id
    predictions <- TransferData(anchorset = lt.anchors, refdata = rc.integrated$time.anno, dims = 1:50, weight.reduction = "pcaproject")
    seu <- AddMetaData(seu, metadata = predictions)
    seu@meta.data$time.anno <- seu@meta.data$predicted.id
    saveRDS(seu, file = paste0("../COPILOT_RDS_LENIENT_DR/",sample.name,"_COPILOT.rds"))
}


#sample.list <- c("sc_12", "sc_17", "sc_35", "sc_57", "sc_58", "sc_59", "sc_60", "sc_61", "sc_62", "sc_63", "sc_64", "sc_66")
sample.list <- c("sc_192", "sc_193", "sc_194", "sc_195", "sc_199", "sc_200", "sc_201", "sc_202")


for (i in 1:length(sample.list)){
 label_transfer_init(sample.list[i])
}

	Meristem	Elongation	Maturation1	Maturation2
	<dbl>	<dbl>	<dbl>	<dbl>
LOC-Os01g01010	3.8171504	3.6033069	3.9269744	4.070318
LOC-Os01g01019	-0.1014923	0.5701927	0.8801323	0.895213
LOC-Os01g01030	3.8802457	3.0066841	3.3899125	3.328385
LOC-Os01g01040	5.3393929	7.0232611	7.1812127	6.560164
LOC-Os01g01050	2.6873846	2.4492419	4.2526756	4.067421
LOC-Os01g01060	9.3628724	7.8553273	6.4893337	5.912760

	Meristem	Elongation	Maturation1	Maturation2	AAACCCAAGAATACAC_1	AAACCCAAGAGTTCGG_1	AAACCCAAGGTAGCAC_1	AAACCCAAGTGATGGC_1	AAACCCAAGTGGAAGA_1	AAACCCAAGTTGCTCA_1	...	TTTGGAGGTTCACGAT_12	TTTGGAGTCCCGTTCA_12	TTTGGAGTCGAACGGA_12	TTTGGTTGTACCTAGT_12	TTTGGTTGTGACTCTA_12	TTTGGTTGTTCGGTCG_12	TTTGGTTTCACCCTTG_12	TTTGTTGAGAAGGATG_12	TTTGTTGAGGATTTAG_12	TTTGTTGCATTCTCCG_12
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
LOC-Os01g01010	3.8171504	3.6033069	3.9269744	4.070318	-0.16539259	-0.1566606	-0.1544584	-0.2539519	-0.2798065	-0.1795019	...	0.36658185	-0.01899874	-0.15512018	-0.20922695	-0.174582428	-0.3825783	-0.1184442	0.11194912	-0.1886373	-0.115706682
LOC-Os01g01019	-0.1014923	0.5701927	0.8801323	0.895213	-0.11405769	-0.1079740	-0.1064359	-0.1748970	-0.1925809	-0.1238417	...	-0.06274370	-0.09550164	-0.04248936	-0.15843951	8.600436939	-0.3794739	-0.1356814	-0.51058476	-0.7206320	-0.072010923
LOC-Os01g01030	3.8802457	3.0066841	3.3899125	3.328385	-0.15887535	-0.1506049	-0.1485180	-0.2424332	-0.2667473	-0.1722251	...	-0.14971078	-0.13922222	-0.13570241	4.30622602	-0.003293015	-0.1495166	5.7988991	-0.32756596	-0.2510060	-0.007144560
LOC-Os01g01040	5.3393929	7.0232611	7.1812127	6.560164	-0.44279801	-0.4109265	-0.4027968	-0.1596220	0.4495326	0.7937728	...	-0.03305562	2.30353379	-0.30791683	-0.04367318	-0.119181227	-0.3277435	-0.1265569	0.03827131	1.5037052	0.007877556
LOC-Os01g01050	2.6873846	2.4492419	4.2526756	4.067421	-0.17543230	-0.1678215	-0.1658997	-0.2519516	2.3777771	-0.1877001	...	-0.18220343	0.14917134	0.06987034	0.94042339	-0.284121127	-0.2629619	-0.2107251	-0.26041142	-0.1918925	-0.129593260
LOC-Os01g01060	9.3628724	7.8553273	6.4893337	5.912760	-0.03083863	-1.1268763	-1.1174055	-1.4572282	-1.5154831	-1.2204200	...	-0.59577470	-0.75737336	-1.04660093	-1.17044143	-0.663788870	-1.4995053	-1.2563877	-0.94665575	-1.1891626	-0.814809660

	AAACCCAAGAATACAC_1	AAACCCAAGAGTTCGG_1	AAACCCAAGGTAGCAC_1	AAACCCAAGTGATGGC_1	AAACCCAAGTGGAAGA_1	AAACCCAAGTTGCTCA_1	AAACCCACACAAGTGG_1	AAACCCACACAGGATG_1	AAACCCACATCAGCGC_1	AAACCCACATTGGCAT_1	...	TTTGGAGGTTCACGAT_12	TTTGGAGTCCCGTTCA_12	TTTGGAGTCGAACGGA_12	TTTGGTTGTACCTAGT_12	TTTGGTTGTGACTCTA_12	TTTGGTTGTTCGGTCG_12	TTTGGTTTCACCCTTG_12	TTTGTTGAGAAGGATG_12	TTTGTTGAGGATTTAG_12	TTTGTTGCATTCTCCG_12
LOC-Os04g46810	-3.489592	-4.076481	-4.064658	-4.107090	-4.281857	-4.181265	-4.3819844	-3.979456	-1.590159	-4.0896639	...	-4.203015	-3.38290498	-1.417949	10.000000	-2.8710460	10.000000	-3.826596	-0.8019649	10.0000000	-2.539364
LOC-Os07g31610	-1.532924	-1.510818	-1.505259	-1.756248	-1.818377	10.000000	-1.7795389	-1.655355	-1.776116	-1.7482597	...	-1.671509	-0.87590624	-1.168881	-1.489915	-1.4419265	-1.770225	-1.648871	-1.8791736	-1.5213221	-1.395994
LOC-Os10g39890	-1.496355	-1.474577	-1.469100	-1.717074	-1.778787	10.000000	-1.7401902	-1.617175	-1.736791	-1.7091515	...	-1.633234	-1.42127324	-1.334040	-1.572713	-1.3881642	-1.689367	-1.650550	-1.5620401	-1.0153073	-1.391670
LOC-Os07g35860	-1.720539	-1.695625	-1.689341	-1.963804	-2.028037	10.000000	-1.9880863	-1.856050	-1.984532	-1.9554223	...	-1.610226	-0.08878541	-1.512589	-1.657318	-1.6891500	-1.895596	-2.163233	-1.8362769	0.2664263	-1.757190
LOC-Os06g22919	-1.510187	-1.487253	-1.481480	-1.739866	-1.802894	10.000000	-0.7884869	-1.636628	-1.760070	-0.6200898	...	-1.713795	-1.23395939	-1.555474	-1.555924	-0.9784831	-1.692611	-1.540442	-1.5317131	-0.7866092	-1.305740
LOC-Os02g44310	-5.422435	-3.098616	-4.656418	-5.437539	-4.585308	-5.607357	-5.6665138	-5.490333	7.762616	-4.6710595	...	-5.739790	-5.38374304	7.562930	10.000000	-3.4914344	10.000000	-3.704694	10.0000000	10.0000000	-1.933807

	orig.ident	nCount_RNA	nFeature_RNA	nCount_spliced_RNA	nFeature_spliced_RNA	nCount_unspliced_RNA	nFeature_unspliced_RNA	percent.mt	percent.cp	nCount_SCT	...	score.anno	consensus.anno	new_ref	final.anno.P	final.anno.cor.P	final.anno.pvalue.P	celltype.anno	time.anno	celltype.fine.anno	integrated_snn_res.500
	<chr>	<dbl>	<int>	<dbl>	<int>	<dbl>	<int>	<dbl>	<dbl>	<dbl>	...	<fct>	<fct>	<chr>	<chr>	<dbl>	<dbl>	<fct>	<chr>	<fct>	<fct>
AAACCCAAGAATACAC_1	tz2	7533	2277	7235	2203	298	141	0.46462233	0.013274924	9654	...	Endodermis	Endodermis	Maturation_Endodermis	Maturation_Endodermis	0.9310545	1.228848e-88	Endodermis	Maturation	Endodermis	341
AAACCCAAGAGTTCGG_1	tz2	6436	1952	6069	1868	367	131	0.52827843	0.000000000	9882	...	Sclerenchyma	Sclerenchyma	Maturation_Sclerenchyma	Maturation_Sclerenchyma	0.8373472	7.925984e-54	Sclerenchyma	Maturation	Sclerenchyma	6
AAACCCAAGGTAGCAC_1	tz2	6178	2392	5902	2268	276	182	0.72839107	0.016186468	9468	...	Exodermis	Cortex	Maturation_Cortex	Elongation_Exodermis	0.4728577	1.548340e-12	Exodermis	Elongation	Exodermis	2067
AAACCCAAGTGATGGC_1	tz2	28987	4895	27735	4684	1252	422	0.42087832	0.000000000	10750	...	Exodermis	Exodermis	Elongation_Exodermis	Elongation_Exodermis	0.8639942	6.486327e-61	Exodermis	Elongation	Exodermis	750
AAACCCAAGTGGAAGA_1	tz2	41734	6475	39910	6181	1824	612	0.06469545	0.002396128	10108	...	Exodermis	Exodermis	Maturation_Exodermis	Maturation_Exodermis	0.8892509	3.514207e-69	Exodermis	Maturation	Exodermis	2602
AAACCCAAGTTGCTCA_1	tz2	9579	3099	9049	2918	530	290	0.48021714	0.010439503	9931	...	Trichoblast	Trichoblast	Maturation_Trichoblast	Maturation_Trichoblast	0.9212941	3.704953e-83	Trichoblast	Maturation	Trichoblast	350

Gene.ID	Cell.type	Source	Atlas.Expression
<chr>	<chr>	<chr>	<chr>
LOC-Os01g14650	Atrichoblast	Zhang et al., Nature Comm., 2021	Confirmed
LOC-Os01g50820	Atrichoblast	Benfey Lab	Confirmed
LOC-Os01g64840	Atrichoblast	Benfey Lab	Confirmed
LOC-Os01g19220	Cortex	Benfey Lab	Confirmed
LOC-Os03g04310	Cortex	Benfey Lab	Confirmed
LOC-Os05g33080	Cortex	Benfey Lab	Confirmed
LOC-Os06g30730	Cortex	Benfey Lab	Confirmed
LOC-Os04g46810	Cortex	Liu et al., Molecular Plant, 2021, Fig. 1E	Confirmed
LOC-Os01g15810	Endodermis	Benfey Lab	Confirmed
LOC-Os01g67390	Endodermis	Benfey Lab	Confirmed
LOC-Os03g18640	Endodermis	Benfey Lab	Confirmed
LOC-Os03g24930	Endodermis	Benfey Lab	Confirmed
LOC-Os10g06680	Endodermis	Benfey Lab	Confirmed
LOC-Os03g02460	Exodermis	Benfey Lab	Confirmed
LOC-Os03g37411	Exodermis	Benfey Lab	Confirmed
LOC-Os04g37980	Exodermis	Benfey Lab	Confirmed
LOC-Os01g16650	Meristem	Zhang et al., Nature Comm., 2021, Fig. S10	Confirmed
LOC-Os02g56130	Meristem	Zhang et al., Nature Comm., 2021	Confirmed
LOC-Os03g51970	Meristem	Zhang et al., Nature Comm., 2021, Fig. S4B	Confirmed
LOC-Os04g41900	Meristem	Zhang et al., Nature Comm., 2021, Fig. S10	Confirmed
LOC-Os05g44400	Meristem	Zhang et al., Nature Comm., 2021, Fig. S4B	Confirmed
LOC-Os01g73980	Metaxylem	Liu et al., Molecular Plant, 2021, Fig. 1E	Confirmed
LOC-Os07g44450	Metaxylem	Zhang et al., Nature Comm., 2021, Fig. S10	Confirmed
LOC-Os01g19170	Pericycle	Benfey Lab	Confirmed
LOC-Os01g58910	Pericycle	Benfey Lab	Confirmed
LOC-Os07g44060	Pericycle	Benfey Lab	Confirmed
LOC-Os07g44060	Pericycle	Zhang et al., Nature Comm., 2021, Fig. S10	Confirmed
LOC-Os04g41350	Phloem	Benfey Lab	Confirmed
LOC-Os06g45410	Phloem	Benfey Lab	Confirmed
LOC-Os08g04400	Phloem	Benfey Lab	Confirmed
LOC-Os08g38170	Protoxylem	Zhang et al., Nature Comm., 2021	Confirmed
LOC-Os03g18130	Root Cap	Ganesh, Alagarasan, et al. 2022 (a review)	Confirmed
LOC-Os03g14300	Root Cap junction	Zhang et al., Nature Comm., 2021	Confirmed
LOC-Os08g02300	Sclerenchyma	Benfey Lab	Confirmed
LOC-Os05g46610	Sclerenchyma	Benfey Lab	Confirmed
LOC-Os01g03530	Trichoblast	Benfey Lab	Confirmed
LOC-Os01g11750	Trichoblast	Benfey Lab	Confirmed
LOC-Os03g04210	Trichoblast	Benfey Lab	Confirmed
LOC-Os04g45290	Trichoblast	Benfey Lab	Confirmed
LOC-Os05g38770	Trichoblast	Benfey Lab	Confirmed
LOC-Os06g08500	Trichoblast	Ding et al., Cell Research, 2009	Confirmed
LOC-Os06g48050	Trichoblast	Benfey Lab	Confirmed
LOC-Os10g03400	Trichoblast	Benfey Lab	Confirmed
LOC-Os10g42750	Trichoblast	Zhang et al., Nature Comm., 2021	Confirmed
LOC-Os01g48130	Xylem	Benfey Lab	Confirmed
LOC-Os01g54620	Xylem	Benfey Lab	Confirmed
LOC-Os01g67090	Xylem	Benfey Lab	Confirmed
LOC-Os09g25490	Xylem	Benfey Lab	Confirmed
LOC-Os10g32980	Xylem	Benfey Lab	Confirmed

Rice Integration¶

Integration¶

Run PCA, UMAP and Clustering¶

Plotting¶

Time zone anno¶

Sample correlation¶

Load the integrated objects (2023/11/09)¶

Cell type anno with known marker z-score¶

Consensus Annotation¶

Correlation annotation round1¶

Correlation annotation round 2¶

Compare cell type /time zone proportion¶

CytoTRACE¶

3D UMAP¶

Save¶

Label Transfer¶